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Abstract. The calculation of distances is of fundamen- 
tal importance in extragalactic astronomy and cosmology. 
However, no practical implementation for the general case 
has previously been available. We derive a second-order 
differential equation for the angular size distance valid 
not only in all homogeneous Friedmann-Lemaitre cosmo- 
logical models, parametrised by Ao and Qq, but also in 
inhomogeneous 'on-average' Friedmann-Lemaitre models, 
where the inhomogeneity is given by the (in the general 
case redshift-dependent) parameter r\. Since most other 
cosmological distances can be obtained trivially from the 
angular size distance, and since the differential equation 
can be efficiently solved numerically, this offers for the first 
time a practical method for calculating distances in a large 
class of cosmological models. We also briefly discuss our 
numerical implementation, which is publicly available. 

Key words: cosmology: theory - methods: numerical - 
cosmology: distance scale - cosmology: miscellaneous - 
gravitational lensing 



1. Introduction 

The determination of distances is one of the most impor- 
tant problems in extragalactic astronomy and cosmology. 
Distances between two objects X and Y depend on their 
redshifts z x and z y , the Hubble constant Ho, the cosmo- 
logical constant Ao, the density parameter Oo and the in- 
homogeneity parameter rjn Usually, smaller distances are 
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1 When discussing the distance between two objects, one can 
always make a coordinate transformation such that the contri- 
bution from the 9 and <f> terms in Eq. (1) vanish. Then one sim- 
ply needs the redshifts and cosmological parameters in order 



determined by the traditional 'distance ladder' technique 
and larger distances are calculated from the redshift, as- 
suming some cosmological model. Since the redshift is for 
most purposes exactly measurable, knowledge of or as- 
sumptions about two of the factors (a) Hubble constant, 
(b) other cosmological parameters and (c) 'astronomical 
distance' (i.e. ultimately tied in to the local distance scale) 
determines the third. In this paper we discuss distances 
given the Hubble constant Hq, the redshifts z x and z y and 
the cosmological parameters Ao, Oo and ij. Traditionally, a 
simple cosmological model is often assumed for ease of cal- 
culation, although the distances thus obtained, and results 
which depend on them, might be false if the assumed cos- 
mological model does not appropriately describe our uni- 
verse. A general method allows one to look at cosmological 
models whether or not they are easy-to-calculate special 
cases and offers the possibility of determining cosmologi- 
cal distances which are important for other astrophysical 
topics once the correct cosmological model is known. 

We stress the fact that the inhomogeneity can be as 
important as the other cosmological parameters, both in 
the field of more traditional cosmology and in the case of 
gravitational lensing, where, e.g. in the case of the time 
delay between the different images of a multiply imaged 
source, the inhomogeneity cannot be neglected in a thor- 
ough analysis (Kayser & Refsdal 1983 ). For an example 
involving a more traditional cosmological test, Pcrlmutter 
et al. ( [1995D (see also Goobar & Perlmutter ( |1995| )) dis- 
cuss using supernovae with z k, 0.25-0.5 to determine go; 
for z near the top of this range or larger, the uncertainty 
due to our ignorance of rj is comparable with the other 
uncertainties of the method. 



to determine the distance between them. When discussing the 
distances between several objects, for example QSOs with a, S 
and z as coordinates, this is no longer possible. In many cases, 
however, suitable geometrical approximations can be made so 
that the most complicated part of the problem is essentially a 
determination of a distance between two objects. This point is 
further discussed in Sect. |5[ 



2 



Kayser et al.: Cosmological distances 



The plan of this paper is as follows. In Sect. ^| the 
basics of Friedmann-Lemaitre cosmology are briefly dis- 
cussed; this also serves to define our terms, which is im- 
portant since various conflicting notational schemes are in 
use. (For a more thorough discussion using a similar nota- 
tion see, e.g., Feige (1992).) Section p[defmes the various 



distances used in cosmology. In SectT|I] our new differen- 
tial equation is derived. Similar efforts in the literature are 
briefly discussed. Section || briefly describes our numerical 
implementation and gives the details on how to obtain the 
source code for use as a 'black box' (which however can be 
opened) for use in cosmology and extragalactic astronomy. 
The symmetry properties of the angular size distance, an- 
alytic solutions and methods of calculating the volume 
element are addressed in three appendices. 

2. Basic theory 

Considering for the moment homogeneous Friedmann- 
Lemaitre cosmological models, we can write the familiar 
Robertson- Walker line element: 



ds 2 = c 2 dt 2 -R 2 (t) 
da 2 



a 2 d0 2 



a 2 sin 2 



where the symbols are defined as follows (with the corre- 
sponding units): 

s 4-dimensional interval [length] 

c speed of light [velocity] 

t time [time] 

R scale factor [length] 

a radial coordinate [dimcnsionless] 

k curvature constant [dimensionless] 

9 angular coordinate [dimensionless] 

4> angular coordinate [dimensionless] 

The dynamics of the universe is given by the Friedmann 
equations 



8nGp{t)R 2 {t) AR 2 {t) 



R 2 (t) = 
and 

R{t) _ AnGp{t) 



- kc 2 



A 

3' 



(2) 



(3) 



R(t) 3 

where dots denote derivatives with respect to t, G is the 
gravitational constant, p(t) the matter density (this pa- 
per assumes negligible pressure), A the cosmological con- 
stant and the sign of k determines the curvature of the 
3-dimensional space. 

Introducing the usual parameters 



* = ! 

R 



n 



A = 



8ttG P 

3H 2 
A 



3H 2 



(Hubble parameter) 
(density parameter) (4) 
(normalised cosmological constant) 



(fl and A are dimensionless and H has the dimension t ) 
we can use Eq. (||) to calculate 



kc 2 = R 2 H 2 (fi + A - 1) , 
so that 

k = sign (O + A - 1) . 
Since R > we can write 

r 1 

R = 



Hy/\n + \-i\' 



(5) 



(6) 



(7) 



this is the radius of curvature of the 3-dimcnsional space 
at time t. For k — it is convenient to define the scale 
factor R to be c/H. In the following the index will be 
used to denote the present value of a given quantity, fixed, 
as usual, at the time to of observation.! 2 ] The explicit de- 
pendence on t will be dropped for brevity. Taking matter 
conservation into account and using the present-day val- 
ues, we have 

pR 3 = p R 3 (8) 
and so from Eqs. @, (4), @ and (||) follows 



/. - H Q R [ — — 



(Qq + A - 1) 



(9) 



Since below we want to discuss distances as functions of 
the cosmological redshift z, by making use of the facts that 



and that Rq is fixed, we can use Eq. (^|) to get 
dz ■ 



dz = — Rdt = -H (l + z)y/Q(z)dt, 

Cl-fl 



(10) 



(11) 



where 

Q{z) = n (l + zf - (r> + A - 1)(1 + zf + A . (12) 

Note: Throughout this paper, the sign should 

be taken to signify the positive solution, except that 
sign y/Q{z) = sign(i?) always. 

3. Distance measures 

3.1. Distances defined by measurement 

In a static Euclidean space, one can define a variety of 
distances according to the method of measurement, which 
are all equivalent. 



2 Note that this paper is concerned with the calculation of 
distances from redshift. We are not concerned with a change 
in redshift with to- 
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3.1.1. Angular size distance 



3.2. Cosmological distances 



Let us consider at position y two light rays intersecting at 
x with angle 9. If I is the distance between these light rays, 
it is meaningful to define the angular size distance D xy as 



D 



xy 



(13) 



since an object of projected length I at position y will 
subtend an angle 9 = l/D xy (for small 9) at distance D xy . 

3.1.2. Proper motion distance 

The proper motion distance is similar to the angular size 
distance, except that I is given by vt, where v is the tan- 
gential velocity of an object and t the time during which 
the proper motion is measured. 

3.1.3. Parallax distance 

Parallax distance is similar to the proper motion distance, 
except that the angle 7r is at y instead of a;, so that we 
have 

D\y = (14) 

In the canonical case, I = 1 AU. 



3.1.4. Luminosity distance 

Since the apparent luminosity L of an object at distance 
D is proportional to 1/D 2 , one can define the luminosity 
distance as 



(15) 



where Lq is the luminosity at some fiducial distance Dq . 



3.1.5. Proper distance 

By proper distance D p we mean the distance measured 
with a rigid ruler. 



3.1.6. Distance by light travel time 

Finally, from the time required for light to traverse a cer- 
tain distance, one can define a distance D c by 



D c = ct 

where t is the so-called look-back time. 



(16) 



3.2.1. General considerations 

In a static Euclidean space, which was used above when 
defining the distances through a measurement descrip- 
tion, these distance measures are of course equivalent. In 
the general case in cosmology, where the 3-dimensional 
space need not be flat (fc = 0) but can be either positively 
(k = +1) or negatively (k = —1) curved, and where the 
3-dimensional space is scaled by R(t), not only do the dis- 
tances defined above differ, but also (in the general case) 
D X y 7^ D yx . The definitions are still applicable, but dif- 
ferent definitions will result in different distances. 

In reality, of course, the universe is neither perfectly 
homogeneous nor perfectly isotropic, as one assumes when 
deriving Eq. (1). However, as far as the usefulness of the 
Friedmann equations in determining the global dynamics 
is concerned, this appears to be a good approximation. 



(See, for example, Longair (1993) and references therein 



for an interesting discussion.) The approximation is cer- 
tainly too crude when using the cosmological model to 
determine distances as a function of redshift, since the 
angles involved in such cases can have a scale comparable 
to that of the inhomogeneities. In this paper, we assume 
that these inhomogeneities can be sufficiently accurately 
described by the parameter 77, which gives the fraction of 
homogeneously distributed matter. The rest (1 — 77) of 
the matter is distributed dumpily, where the scale of the 
dumpiness is by definition of the same order of magnitude 
as the angles involved. 

For example, a halo of compact MACHO type objects 
around a galaxy in a distant cluster would be counted 
among the homogeneously distributed matter if one were 
concerned with the angular size distance to background 
galaxies further away, but would be considered clumped 
on scales such as those important when considering mi- 
crolensing by the compact objects themselves. Since we 
don't know exactly how dark matter is distributed, dif- 
ferent 77 values can be examined to get an idea as to how 
this uncertainty affects whatever it is one is interested in. 
If one has no selection effects, then, due to flux conser- 
vation, the 'average' distance cannot change (Weinberg 



1976); 77 introduces an additional uncertainty when inter- 



preting observations. It is generally not possible to esti- 
mate this scatter by comparing the cases 77 = and 77 = 1, 
since, depending on the cosmological parameters and the 
cosmological mass distribution, not all combinations are 
self-consistent. For instance, if one looks at scales where 
galaxies are compact objects, and the fraction of f2o due 
to the galaxies is x, then 77 must be < (1 — x). 

We further assume that light rays from the object 
whose distance is to be determined prop agate sufficiently 
far from all clumps. (See Schneider et al. ( 1992] ) - hereafter 
SEF - for a more thorough discussion of this point.) Com- 
pared to the perfectly homogeneous and isotropic case, the 
introduction of the 77 parameter will influence the angular 
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size and luminosity distances (as well as the proper mo- 
tion and parallax distances) since these depend on angles 
between light rays which are influenced by the amount of 
matter in the beam, but not the proper distance and only 
negligibly the light travel time. The last two distances are 



discussed briefly in Sect. 3.2.2 and in App. B.3 and B.6 



Since there is a simple relation between the angular size 
distance and the luminosity distance (Sect. 3.2.2] ) which 
also holds for the inhomogeneous case (see App. |a]), for 
the general case it suffices to discuss the angular size dis- 
tance, which we do in Sect. [|. 

3.2.2. Relationships between different distances 

Without derivation! we now discuss some important dis- 
tance measures, denoting the redshifts of the objects with 
the indices x and y. Due to symmetry considerations (sec 
App.0) 



D„ 



l + z v 
1 + z x 



(17) 



where the term in parentheses takes account of, by way of 
Eq. ([lOj), the expansion of the universe. It is convenient, 
in keeping with the meaning of angular size distance, to 
think of the expansion of the universe changing the angle 
in Eq. ( |l3| ) and not I, if one identifies I as the (projected) 
size of an object. The angle is defined at the time when 
the light rays intersect the plane of the observer. Thus 
D xy with the observer at x = defines what one normally 
thinks of as an angular size distance. On the other hand, 
D xy and D yx with x in general 7^ can be important in, 



for example, gravitational lensing.l 

Although the angle between the rays (at the source) 
at the time of reception of the light is important for the 
luminosity distance, this distance is not simply D yx , since 
in the cosmological case the observed flux is obtained 
by multiplying the 'non-redshifted flux' by the factor 
(l+z x ) 2 /(l+z y ) 2 . One factor of (l+z x )/(l+z y ) occurs be- 
cause a given wavelength is increased by (1 + z y )/{ \ + z x ), 
which reduces the flux correspondingly; an additional fac- 
tor of (1 + z x )/(l + z y ) occurs because the arrival rate 
of photons is also decreased. Therefore, since D L is in- 
versely proportional to the square root of the (observed, 
'redshifted') flux the luminosity distance is 



(18) 



3 See, e.g., Feige ( |l992| ) Berry ( |1986[ ) or Bondi (g_96l|) for a 
more general discussion. What we present in the rest of this 
section is not new, but is important in order to clarify the no- 
tation. The results are obvious from the definitions introduced 
above. 

4 Although not useful in cosmology or extragalactic astron- 
omy, for completeness we mention the fact that the proper 
motion distance is equivalent to D yx and the parallax distance 
is equivalent to D yx /yjl — ka 2 . 



From this and Eq. (17) follows the relation 
1 + ; 



^ = D X y 



1 



(19) 



This means that the surface brightness of a 'standard can- 
dle' is ~ (1 + z)~ 4 , a result independent of the cosmologi- 
cal model parameters, including r/B (This result also holds 
for the inhomogeneous case, since Eq. (17) still holds (see 
App. [X]) and the additional factor due to the expansion of 
the universe (given by the term in parentheses in Eq. (|L8|)) 
is of course present in the inhomogeneous case as well.) 

Of course, this applies only to the bolometric lumi- 
nosity. Observing in a finite band introduces two correc- 
tions. The so-called ^-correction as it is usually defined 



today (see, e.g., Coleman et al. (1980) or, for an interesting 
and thorough discussion, Sandagc (1995)) takes account 
of these, both of which come from the fact that the ob- 
served wavelength interval is redshifted compared to the 
corresponding interval on emission. This means that, first, 
for a flat spectrum, less radiation is observed, because the 
bandwidth at the observer is (1 + z) times larger than 
at the source. Second, the spectrum need not be flat, in 
which case additional corrections phased on the shape of 
the spectrum have to be included!!] Thus, 



m = M + 5 log 



10 pc 



K 



(20) 



where m is the apparent magnitude, M the absolute mag- 
nitude, D L is the luminosity distance and K is the K- 
correction as defined in Coleman et al. (1980). Perhaps 
more convenient is 



M + 5logD h +K + N 



(21) 



where N is a normalisation term: N = — 5 for D L in units 
of 1 pc, N = 25 for D L in units of 1 Mpc jpd N = x — 
51og/i for D L in units of the Hubble lengthQ c/Hq, where 



x 



5 log 



/ Hubble length 



1 pc 



5 « 42.384 



5 Thus, a 'surface brightness test' can in principle show that 
cosmological redshifts are due to the expansion of the uni- 
verse and not to some other cause. See, e.g., Sect. 6 in Sandage 
(jl995| ). 

Since the observed objects generally evolve with time, and 
redshifted objects are necessarily observed as they were when 
the radiation was emitted, some authors include an evolution- 
ary term in the Tf-correction. Still other authors prefer to 
absorb one or more of these terms into the definition of the 
luminosity distance. Our luminosity distance is a bolometric 
distance based on the geometry and includes the unavoidable 
dimming due to the redshift. Our A'-correction takes account 
of both effects of a finite bandwidth. Evolutionary effects are 
considered separately from distances. 

7 For example, as given by our numerical implementation; 
see Sect. |H| 
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and h is the Hubble constant in units of lOOkm/s/Mpc. In 
practice one has to add terms to correct for various sources 
of extinction and consider the fact that M is the absolute 
magnitude of the object when the light was emitted, which 
of course could be different from the present M of similar 
objects at negligible redshift. 

The light travel time (or lookback time) t xy = t x — t y 
between z x and z y (where t x = t(z x ) > t y — t(z y )) is given 
by the integration of the reciprocal of Eq. (0) : 



dz 



dz 



H J (l + z)^Q(i)' 



(22) 



two light rays on the boundary of a small light cone prop- 
agating far away from all clumps of matter in an inhomo- 
geneous universe: 



R 

I = -AnG-qp 1 + —l 
R 



(24) 



where r\ and p are functions of the time t (not the lookback 
time of Eq. ^2|). The first term can be interpreted as Ricci 
focusing due to the matter inside the light cone, and the 
second term is due to the expansion of space during the 
light propagation. We now have to transform this time 
dependent differential equation into a reds/ii/i-dependent 
differential equation. From Eq. (O) we obtain^ 



where the minus sign from Eq. is equivalent to the dt = — \Hq(1 + z)\/Q) dz 
swapped limits of integration on the right-hand side so 
that the integral gives t x — t y instead of t y — t x , making 
the light travel time increase (for R > 0) with z; thus 



and thus 



D 



Ct^f 



Since the proper distance would be the same as D c 
were there no expansion, the former can be calculated by 
multiplying the integrand in Eq. (|22|) by c(l + z). Thus 



= — 



dz 



(23) 



This gives the proper distance at the present time. Since 
D F scales linearly with the expansion of the universe, the 
proper distance at some other time can be obtained by di- 
viding Eq. ( p3| ) with (1+Zi), where Zi is the redshift at the 
corresponding time. For homogeneous (rj = 1) cosmolog- 
ical models,u the propagation of light rays is determined 
by the global geometry, so that there is a simple relation 
between D v and D and, thus, D h . This is discussed in 



Sect. [B.3| . Although not 'directly' observable, the proper 
distance is nevertheless important in cosmological theory, 
since it is the basic distance of general relativity. Although 
not useful as a distance, the light travel time is of course 
important when considering evolutionary effects. 

For inhomogeneous models, where this relation be- 
tween global geometry and local light propagation does 
not exist, another approach must be used, which takes ac- 
count of both the expansion of the universe as well as the 
local propagation of light, when calculating angle-defined 
distances such as the angular size distance. 

4. The general differential equation for the angular 
size distance 



In a seri es of papers Zeldovich (1964), Da shevs kii and Zel- 
dovich ( 1965|) and Dashevskii and Slysh ( 1966 ) developed 
a general differential equation for the distance between 



8 This includes empty models (fio — 0); although r\ has no 
meaning here, the same arguments apply. 



dl tt h n m. dl 

and 



(25) 



(26) 



i*(((i + „g + (i + ,>-ig)£ 



dz z 



(28) 



Furthermore, since R = Rq/(1 + z) (Eq. fllCD), we obtain, 
using Eq. ©, 



(29) 



— = -H (l + z)^Q — 
dt dz 



From the definition of SI (Eq. (4)) and matter conservation 
(Eq. d§)) we obtain 



4ttG P = -Hfooil + zf 



(30) 



If we now insert Eqs. (§|), (28), © and Q into 
Eq. ([24]), sort the terms appropriately and cancel Hq, 
which appears in all terms, we obtain 

Ql " + {tT~ z + I q/ ) l ' + lvW + z)i = o , (31) 

where a prime denotes a derivative with respect to redshift 
and from Eq. (Il2]) follows 



(32) 



Q'(z) = 3fl (l + z? ~ 2(ft + A - 1)(1 + z) 



9 This transformation causes problems if the integration in- 
terval contains a point where R — and thus \[Q changes 
sign. In this case the integration interval (t x ,t y ) has to be 
transformed into two integration intervals, namely (z x , z max ) 
and (z max ,z y ), where z max is the redshift at R — 0, with the 
boundary conditions for the second integration interval chosen 
appropriately. 



(i 
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From the definition of the angular size distance (Eq. (|13[)) 
it is obvious that it follows the same differential equation 
as Z: 



QD" 



2Q 

1 + 2 



1 



Q'\ D' + - 7 1 fl (l + z)D = 



(33) 



with special boundary conditions at the redshift z x where 
the two considered light rays intersect. The first boundary 
condition is trivially 



D = for 



(34) 



and the second boundary condition follows from the Eu- 
clidean approximation for small distances, i.e. 



AD 

hence 
D' 



csign^ -t y ), 



(35) 



Ho (1 + Z x )y/Q&) 



siga(t y -t x ) for z = z x , (36) 



where the sign has been chosen such that D is always 
> locally. We denote these special solutions of Eq. ( j33"[) 
with D x (z), and, following the definition (Eq. (|l3|)), the 
angular size distance of an object at redshift z y is then 



given as 



D X y D X {Zy 



(37) 



Figure [l] shows the influence of z, r\ and A on the an- 
gular size distance, calculated using Eq. ( p3| ) with our nu- 
merical implementation. 

For completeness we note that after the original deriva- 
tion by Kayser ( 1985| ) an equivalent equation was derived 
by Linder (1988) which, however, is difficult to implement 
due to the cumbersome notation. 

Special mention must be made of the so-called bounce 
models, which expand from a finite R after having con- 
tracted from R = oo. (See, e.g., Feige ( 1992 ).) A glance 
at Eq. ([Io|) shows that in these cosmological models there 
must be four distances for an (ordered) pair of redshifts. 
If we denote the distances by D\2, -D14, -D34 and -D32, 
where 1(2) und 3(4) refer to 24(2:2) during the expanding 
(contracting) phase, then symmetry considerations dictate 
that D12 = -D34 and D14 = -D32 as long as the dependence 
of Tj on z is the same during both phases. In this case, 
there are two independent distances per (ordered) pair 
of redshifts. If this is not the case, the degeneracy is no 
longer present and there are four independent distances 
per (ordered) pair of redshifts. 

5. Numerics and practical considerations 

For the actual numerical integration of the differential 
equation, we have found the Bulirsch-Stoer method to be 




Fig. 1. The angular size distance from the observer (21 = 0) 
and from 21 = 2 (lower right) as a function of the redshift Z2 
for different cosmological models. Thin curves are for n — 0, 
thick for r\ = 1. The upper curves near 2 = (2 = 2 at lower 
right) are for Ao = 2, the lower for Ao = 0. Qo — 1 for all 
curves. The angular size distance D is given in units of c/Ho 



both faster and more exact than other methods such as 
Runge-Kutta. However, the conventional method of ratio- 
nal function extrapolation is rather unstable in this par- 
ticular case; fortunately, using polynomial extrapolation 
solves the problem. Although programming the integra- 
tion is rather straightforward in theory, in numerical prac- 
tice considerable effort is needed to determine combina- 
tions of free parameters which work for all cases. We have 
tested the finished programme intensively and extensively, 
for example by comparing the results of calculations for 
n = 1 (the value of 77 plays no special role in the integra- 
tion ofthe differential equation) with those in Refsdal et 
al. ( 1967 ) or give n by t he method of elliptical integ rals a s 
outline d in Fe ige ( 1992 ) and have use d it in Kayser ( 1995 ), 
Helbig ( 1996 ) and Helbig & Kayser ( |1996|) . For a general 
discussion of various methods of integrating second-order 
differential equations, see Press et al. ( |1992| ). Those inter- 
ested in technical details can read the comments in our 
source code and the accompanying user's guide. 

Since Ho, in contrast to the other cosmological param- 
eters, merely inversely scales the angular size distance, 
our routine actually calculates the angular size distance 
in units of c/H Q . This dimensionless quantity must be 
multiplied by c/ Hq (in whatever units are convenient) in 
order to obtain the actual distance. Other than reducing 
numerical overhead, this allows all distances to be cal- 
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culatcd modulo c/Hq, which is convenient for expressing 
quantities in an i/o-independent manner. In practice, Hq 
cancels out of many calculations anyway. 

Apart from auxiliary routines which the user does not 
have to be concerned with, our implementation consists of 
four FDRTRAN77 subroutines. The first, INICOS, calculates 
z-independent quantities used by the other routines, some 
of which are returned to the calling programme. ANGSIZ 
calculates the angular size distance. Normally, 77 is used 
as a z-independent cosmological parameter, on an equal 
footing with Ao and fio- If desired, however, the user can 
let INICOS know that a variable (that is, z-dependent) 77 is 
to be used; this is given by the function VARETA. We supply 
an example; the user can modify this to suit her needs. In 
particular, many different dependencies of 77 on z can be 
included, and a decision made in the calling programme 
about which one to use. This feature is also included in 
our example. ANGSIZ returns only the distance D12; if one 
is interested in the other distances in the bounce mod- 
els, our subroutine BNGSIZ returns all of these (though 
internally calculating only the independent distances, of 
course, depending on the dependence of r\ on z). 

Due to the fact that not everyone has a Fortran90 
compiler at his disposal, we have coded the routines in 
F0RTRAN77. Only standard FORTRAN 7 7 features are used, 
and thus the routines should be able to be used on 
all platforms which support F0RTRAN77. Since standard 
F0RTRAN77 is a subset of Fortran90, the routines can be 
used without change in Fortran90 as well. 

With the exception of D c , all distance measures can 
be easily transformed into one another. Thus, it suffices 
to calculate the angular size distance for a given case.Ej 

When discussing the distance between two objects 
other than the observer, rather than between the observer 
and one object, in many cases one of two simplifying as- 
sumptions can be made: 

D(Az) <C D{(3) In this case, the proper distance D p at the 
time of emission between the two objects is (3D$ X « 
[3Do y , where (3 <C 1 is the angle in radians between the 
two objects on the sky. 

D(f3) <C D{Az) In this case, the angular size distance be- 
tween the two objects is D xy . 

D(Az) (D((3)) refers to the distance due to Az ((3) when 
setting (3 (Az) equal to zero. In the first case, where the 
two objects are practically at the same redshift, one uses 
the angular size distance to this redshift to transform the 
observed difference in angular position on the sky into the 
proper distance between the two objects at the time of 
emission. This follows directly from the definition of the 



angular size distance. Since the distance between the ob- 
jects is much less than the distance from the observer to 
the objects, the differently defined distances between the 
objects are for practical purposes degenerate. A practical 
example of this case would be the distance between indi- 
vidual galaxies in a galaxy cluster at large redshift. Natu- 
rally, one should use one redshift, say, of the cluster centre; 
the individual redshifts will in most cases be overlaid with 
the doppler redshift due to the velocity dispersion of the 
cluster, so the difference in cosmological redshifts is neg- 
ligible. (Of course, the present distance would be a factor 
of (1 + z) larger, due to the expansion of the universe, 
were the objects comoving and not, as in a galaxy cluster, 
bound.) In the second case, which is typical of gravita- 
tional lensing, the angles on the sky between, for exam- 
ple, source and lens, are small enough to be neglected, so 
that the angular size distance between the objects is de- 
termined by the difference in redshift. If neither of these 
assumptions can be made, any sort of distance between 
the two objects is probably of no practical interest. (Of 
course, there is the trivial case where the redshifts are all 
<C 1 in which case one can simply use a, S and cz/Hq as 
normal spherical coordinates.) 



6. Summary 

After discussing cosmological distances with an emphasis 
on practical distance measures for general use in cosmol- 
ogy and extragalactic astronomy, we have obtained a new 
differential equation, which gives the angular size distance 
for a class of 'on average' Friedmann-Lemaitre cosmolog- 
ical models, that is, models described not only by Ao and 
f^o but also by ?7(z), which describes the dumpiness of the 
distribution of matter. We have also developed a practical 
numerical method of solving this equation, which we have 
made publicly available. Since the equation is valid for all 
cases, this offers for the first time an efficient means of cal- 
culating distances in a large class of cosmological models. 

The numerical implementation (in F0RTRAN77), user's 
guide and a copy of the latest version of this paper can be 
obtained from either of the following URLs: 



tittp: / /www. hs.uni-hamburg.de/cnglisli/persoiis/hclbig/ 
Research / Pu blications /Into/ angsiz. html 



Ftp: / / ftp.uiii-hamburg.de/pub/uiiihh/astro/aiigsiz.tar.gz 



10 The proper distance, which is ^-independent, can be calcu- 
lated from the angular size distance assuming n — 1, by making 
use of the simple relation between proper distance and angu- 
lar size distance in this case. The result holds of course for all 
values of n. 
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A. Symmetry: The relation between D xy and D yx 

The proof in this appendix follows closely the proof pre- 
sented in Kayser ( 1985 ). For completen ess we note that 
after the original derivation by Kayser (1985) an equiv- 
alent equation was derived by Linder (1988). We rewrite 
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the differential equation, Eq. 
tance in the normal form: 



for the angular size dis- and inserting a, and fa back into Eq. ( jAq ) we obtain for 

the special solutions D x : 



a 2 D"(z) + oi(«) £>'(z) + a Q (z) D(z) 
with the coefficient functions 
a 2 (z) = Q{z) 







a\{z) 



2Q(z) 
l + z 



1 



Q'(z) 



(Al) 

(A2) 
(A3) 

(A4) 



D 



(A14) 



If wc now consider a second special solution D y we find 
the relation 



(A15) 



00(2) = 2 ^^0(1 + 2) ■ 

The Wronskians can be calculated using Liouville's for- 
Now let and L> (2) be two solutions of Eq. jAl| ) which mula: 
build a fundamental system, i.e. the Wronskian for these 
two solutions does not vanish: 



W(z) 



D W £,(2) 



dz 



dz 



(A5) 



Every solution Di of Eq. (Al) can then be written as a 
linear combination of and D^: 



D, 



a>i D (1) + fa L> (2) ; with a i} fa = const 



The angular size distances are special solutions D x of 
Eq. (Al) fulfilling the following boundary conditions: 



D. r = for z = 



and 

dD x 
dz 

with 



b(z x ) for z = z x 



b{z x ) = -Jr- ((1 + z x )^/Q(z x )j sign(t y - t x ) 



compare Eq. (pq). From Eq. (A6) we obtain 

Q = a t D^{z x )+faD^{z x ) 
and 



(A7) 



(A8) 



(A9) 



(A10) 



b(z x 



dz 



+ A 



dZ?( 2 ) 



dz 



These equations can easily be solved for a, and fa: 

D^(z x ) 



Oti = (3, 



DW(z x ) 

b(z x )D^(z x ) 
W{z x ) 



W(z) = W(z ) exp a 2 {z)dz , 



(A16) 



where zq is arbitrary. Thus 



Dx(Zy) 

(A6) D y (z x ) 



b{z x ) 

b(Zy) 



exp 



ai(z) 
a 2 (z) 



dz 



(A17) 



and after inserting ao, a\ and a2 from Eqs. (A2), ( |A3| ) 
and (A4) as well as b(z x ) and b(z y ) from Eq. (A9) and 
integration we finally obtain for the angular size distances 
(cf. Eq. @) the relation 



D, 



1 + z x 

l + Zy 



(A18) 



B. Special cases 

For certain special cases the differential equation can be 
simplified and sometimes analytically solved. 



B.i. a, = 



A glance at Eq. ( |33[) shows that for O = the third 
(All) term on the left hand side of Eq. ( |3^ ) vanishes; one thus 
has a first order differential equation for D' . (Of course 
j] has no meaning for £!o = 0.) Due to the fact that a 
vanishing f^o also simplifies Q(z), it is possible to calcu- 
late the angular size distance analytically. Since in this 

(A12) 

case the angular size distance is determined exclusively 
by global (effects, one can use an approach based on global 
geometry.Eil Depending on the value of Ao, one can use the 

(A13) 



See the discussion in Sect. 



B.3. 
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following expression to calculate Xxy = x(y) — x( x ) (Feige B.3. 77 = 1 



1992) 



xO) = < 



'arccosh(V') for A < 

ln(l + z) for A = 

arcsinh(V') for < Ao < 1 

z for Ao = 1 

arcsin(V>) for Ao > 1 



(Bl) 



where tp := (1 + z )y |Ao|°^ ' relationship between % 
and the angular size distance D is 



sinhx for k = — 1 
X for k = 
+ ^ I sinx for k = +1 



(B2) 



as discussed below in Sect B.3 



r? = 

In the case n = the third term on the left hand side 
of Eq. ( |33|) vanishes; one thus has a first order differential 
equation for D' . Assuming D' ^ 0, Eq. ( |33"1) can be written 

as 



(B3) 



D" _ 2 1 Q'(z) 
~D> ~ ~l + z ~ 2 Q(z) ' 

This equation can be solved in two steps. For D' we obtain 



The case r\ = 1 has all matter distributed homogeneously. 
Due to homogeneity, the matter locally affecting the prop- 
agation of light is known when the global geometry is 
known, so that the 'classical' approach of relating global 
geometry to observable relations is a better approach than 
using (the simplified form of) Eq. ([53]) . This approach of- 
fers an analytic solution. Here, we simply sketch the most 
i mpor tant points; the interested reader can refer to Feige 
( 1992 ) for a good description of this method. 
The angular size distance in this case is 



xy ~ y(Jxy ~ (1 + z v ) ' 



(B9) 



where a is the radial coordinate in the Robertson- Walker 
metric (cf. Eq. (1)) and thus 



Rncr : 



1 + Z X 



= Dxy 



l + Z V 
1 + z x J ' 



(BIO) 



since this angle is inversely proportional to R for constant 
a and physical size. (The value of R at the time the light 
rays defining the angle intersect is important.) 
Since a is given by 



sinhx for k = — 1 
o = n v ! = { X for k = , 
sinx for k = +1 



(Bll) 



D' 



ci 



and consequently for D 
D= ' 



c 2 . 



(B4) 



(B5) 



The constants c\ , C2 are determined by the appropriate 
boundary conditions (Eqs. ( |34| ) and (|3^)). We then find 
the solution (see also SEF for an equivalent discussion 
with Ao = 0) 



D xy = — (1 + z x )(uj(z y ) - uj(z x )), 
Ho 



where 
uj(z) = 



dz' 



(1 + z')V(i + ^) 2 (fV + l-A ) + Ao' 



or, perhaps more convenient, 

c , . f Zy dz 

D x V = jr( 1 + z *) / 77~ xT7nr\ 



(B6) 



(B7) 



(B8) 



For Ao = there is an analytic solution (see Sect. B.4 ) 



an expression for x( z ) is sufficient for calculating the an- 
gular size distance D (and of course the luminosity dis- 
tance D L (via Eq. (pj|)) and the 'coordinate distance' a 
(via Eq. (Bll)). In general, a xy ^ o~ y — a x ; however, 



Xxy — Xy Xx, 



F(Xxy) 



so that 



(B12) 



where F is given by Eq. (Bll). Using Eq. (23) one can 
calculate 



Xxy 



El 

Rq 



HqRq 



dz 



(B13) 



In the general case, Eq. (23) can be solved by elliptic in- 
tegrals, as explained in Feige (1992). For the cases Ao = 
and f2o = the formulae using elliptic integrals break 
down; in these cases, easier analytic formulae, which for- 
tunately exist, can be used. The case JIq = has been 
discussed above. The case A = will be discussed below. 
Again, we stress that the differential equation derived in 
Sect. U is completely general and can be used in all cases. 
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B.l A = 



B.4.2. An = and O n = 1 



For An — 0, there is in general no simpler solution. This 
case h as been discussed by Dyer and Roeder for -q = 
( |1972D and for general r\ values ( 1973 ). They point out 
the interesting result that the maximum in the angular 
size distance from Z\ = to z 2 increases monotonically 
from 1.25 to oo as r\ decreases from 1 to 0. See also the dis- 
cussion (with a differing notation!) in Sect. 4.5.3 in SEF. 
However, some solutions exist for special values of Oo and 
r\. The case fin = has been discussed in Sect. B.l above; 



the value of r\ is of course irrelevant in this case. With the 
exception of Sect B.4.3 below, in the following we simply 
quote results from SEF in our notation. 

B.4.1. A = and r\ = 

As discussed above, for r) = Eq. j33|) is effectively a first 
order equation for D'. For Ao = Q(z) is sufficiently sim- 
plified to allow an analytic solution. Recalling Eq. (B6), 



D, 



H„ 



(1 + Z x )(u(Zy) - U>(z x )), 



Eq. (B7) simplifies to 

; dz' 

uj{z) 







(1 + OVcWTT 



(B14) 



which has the solution 

3fif 



uj{z) 



~ w arctan(^) + AIO , ri 

5f2 n 



|)yo z+i 



(T+zy 2 



1 - {Q z + 



— — — r arctanh(-0) 

4(l-O )3 

(n oZ +Sai-|)vT1^+ 

(1 + z) 2 



3 



(I) 



(II) (B15) 



4(i-n ) 2 



3 



(in) 



with 



and 



v / n^T(i+vn z+i) 
Qo-i+vnoz+i 



VT^n^(i+vQ z+i) 

Oo-l+V^oz+l 



(I) 
(III) 



(B16) 



case I: O > 1 

case II: fin = 1 
case III: < O < 1 



(B17) 



Note that in SEF, the text at the top of page 137 is 
unclear — the expression in parentheses in the denomina- 
tor of the first term (Q — 1) for the O > 1 case has to be 
replaced with (1 — f2) as well for < 1. Note also that 
O = f2o and that after page 131 An = is always assumed. 



For fio = 1 and Ao = (the Einstein-de Sitter model) we 
have the solution 



1 / {l + z y f-i (l + z x f-i 



where 



/3:=-V25 - 24»7. 



B.4.3. A = and 77 = 1 



(B18) 



(B19) 



For r\ = 1 the special case of the expression for x( z ) f° r 
A = is (Fcigc |1992D 



X(z) 



-2 < 



(f! = 1) , (B20) 

(o < < i) 



arcsinh I 



where Xxy = x{y) ~ X. ( x )- (It is obvious that in the case 
A = Ho = Eq. (Bl) should be used.) From this, it is 
possible to obtain a general expression for the angular size 
distance (see, e.g., SEF): 



D *y = 77-7^2 (l+*x) {Ri(z y )R 2 (z x ) - R l {z x )R 2 {z v )) ,(B21) 
n sz 



with 

Ri(z) 
and 

R2(z) 



n z - o + 2 

(1 + z) 2 



yjn a z + 1 

(1 + z) 2 • 



(B22) 



(B23) 



For z x = and z y — z one gets for the angular size dis- 
tance 



H 2 (1 + zf 
(p z - (2 - o ) (y^ 



I lj). (B24) 
valid for Oq > 0. For Oq = one obtains 



D 



c z 



(1+1) 



#0 (i+^r 



(B25) 



(Multiplying Eq. (B24) or Eq. flB25|) with (1 + z)/i? re- 
sults in the respective expression for a as a function of red- 
shift as first derived by Mattig (1958). See also Sandage 
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( |1995| ), Sect. 1.6.3). In this case, the volume element given 
by Eq. (H3) reduces to 



dV = 167Ti? 3 ^° Z " (2 " n 0)(V^+l - I) 2 ) 



ng(i + «)V«o« + i 



(B26) 



Of course, for the physical, as opposed to comoving, den- 
sity, an additional factor of (1 + z) 3 must be added to the 
denominator. 



BAA. A = and r) = § 

For r\ = | and Ao = there is also an analytic solution 
(see SEF): 



*5 S '0 

(R 1 (z x )R 2 {z y ) - R 2 {z x )Ri{z y )) 



with 



Ri(z) = 



and 



R 2 (z) = 



(l + zf 



^/tt z + l[n z + 3O - 2) 



(l + zf 



(B27) 
(B28) 

(B29) 



B.5. Other 



We can offer no proof that no other easier solutions, ei- 
ther reducing Eq. ( |33| ) to a more easily (numerically) in- 
tegrated form or even to an analytic solution, exist. This 
is left as an exercise to the interested reader. The authors 
are of course interested in such solutions and are willing 
to verify them. As far as we know, Eq. ( |33| ) must be used 
except in the special cases mentioned in this appendix. 

B.6. Light travel time 



Feige (|l992j ) not only gives the distance but also the light 
travel time by means of elliptic integrals. As for the dis- 
tance, and for the same reasons, simple analytic formulae 
can and must be used for the special cases Slo = and 
Ao = 0. For fc = 0, an analytic expression for the light 
travel time exists, although the elliptic integrals can also 
be used in this case. For completeness, we give these spe- 
cial cases here for the light travel time t xy = t x — t y . 
For SIq = we have: 



'arcsm(i/>) 



t(z) 



1 



A < 



/ |A„| 



tfoVlAol 



Ao = 
< Ao < 1 . (B30) 



arcsinh(-0) 
-^AoTln(l + z) A = l 
_ arccosh(V>) Aq > 1 



where ip 



[Ao] 

(l+z) V |1— A | 



For An = we have: 



n (i+z) 



arcsm | 



t(z) = A x < 



3 \y i+z 



n > i 

n = i ,(B3i) 



arcsinh 



where 
.4 



- v /(n z+i)(i-n ) 
n (i+2) 



n 



2o(l+z) J 



< fin < 1 



ff (^o^l|) 3 



(For Cl = the appropriate case from Eq. ( B30 ) must be 
used.) 

For k = we have: 



'7n7=T arcsill W 



t(z) 



3H 



■ x < 



T , 1 _ arcsinh(-!/;) 



n Q > i 

Qo = l ,(B32) 
0<il o <1 



where 



n„(i+z)3 



Sl > 1 
< Sl < 1 



(B33) 



(For f^o = the appropriate case from Eq. ( B30 ) must be 
used.) 

C. Volume element 

Sometimes the distance is only a means of calculating the 
volume element at a given redshift. In the static Euclidean 
case the volume element is of course 



dV = 4irr dr. 



(CI) 



In the cosmological case, the volume element is, with r 
R a, 



dV = 4irr 2 dD^ = Airr 2 



(C2) 



RoCTy is, for r) = 1, simply D y0 = (1 + y)D 0y ; see Sect BJ3 



Thus, the distance D y0 is all that is needed to calculate 



the volume; this first can be calculated by Eq. (|33|) with 
rj = 1 (This applies even if one would calculate distances 
with another value of rj since the volume element is a 
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quantity related to the global geometry of the universe- 



alternatively, one can use elliptic integrals, as in Sect. B.3 
and Feige ( 1992f ).) If one has an expression for cr(z), then 
since 



dD p = Rq- 



which follows directly from Eq. (1), one can write 



dV(a) 



where Rq is given by Eq. (^) for the present values: 
R C 1 



Ho v/^o + Ao-ll 



and 



a = D ya( l + Z v) 
Integration gives 



V(a) = { |7rr 3 



o / Vl-\-(T 2 arcsinh a 



3 

2nr 3 



( arcs 

V 17 



csin cr 
■3 



—Jl — 



fc = -1 

fc = 

fc = +1 



(C3) 



(C4) 



(C5) 



(C6) 



(C7) 



-1, the total volume of the universe is 



Thus, for fc 

2vr 2 i?§. (See, e.g., Sandage ( |1995| ), Sect. 1.6.1; Sandage's 
d is our D p and his r is our a.) Since 
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da 



Vl - fcc 



Eq. flC7| ) can also be written as (cf. Feige ( |1992|) 
Eq. (116); Feige's r is our a) 



V(x) = 2nRl { fx 3 



sinh(x) cosh(x) — \ 
2 Y 3 

X - sin(x) cos(x) 



fc = -1 

fc = (C8) 
fc = +1 



Of course, all this refers to volumes now at the distance 
corresponding to z = y. If the volume at another time is 
important, say at the time of emission of the light we see 
now — for instance if one is concerned with the space den- 
sity of some comoving objects — then the volume element 
must be divided by (1 + z) 3 . 
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